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SUMMARY 


An improved two-phase model has been developed for a wiped-fiim (third stage) 
PET reactor. The model accounts for all the important main and side reactions and 
incorporates the effect of vaporization of four low molecular weight volatile species. 
Industrial data under three different operating conditions are used to obtain best-fit 
(tuned) values of ten model parameters. 

Multiple objectives (often conflicting) occur naturally in polymerization reactors 
and offer a great deal of challenge and insight when vector optimization is applied to 
these systems. Polyethylene terephthalate (PET) is one of the most widely used polymers 
in day-to-day life. The profit margins in synthetic fiber industries (especially PET fibers) 
are declining and therefore any economic improvement through process optimization 
would cause a major reduction in the operating costs. In an attempt to achieve this 
objective, we have undertaken a task of fundamental modeling and multiobjective 
optimization of the wiped-film finishing reactor as a preliminary step in the process of 
system optimization. Multiobjective optimization of the third-stage wiped film polyester 
reactor is then carried out using the ‘tuned’ model. The two objective functions 
minimized are the acid and vinyl end group concentrations in the product. These are two 
of the undesirable side products produced in the reactor. The optimization problem 
incorporates an end-point constraint to produce polymer having a desired value of the 
degree of polymerization (DP). In addition, the concentration of the di-ethylene glycol 
end group in the product is constrained to lie within a certain range. Two adaptations of 
the nondominated sorting genetic algorithm (NSGA - I and NSGA - II) have been used 
to solve the above multiobjective optimization problem. The decision variables used are 
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the reactor pressure, temperature, catalyst concentration, residence time of the polymer 
inside the reactor and the speed of the agitator. A unique optimal solution was obtained in 
all cases unlike what was observed in an earlier study [Bhaskar et al. (2000)]. Use of 
suitable GA parameters is necessary to avoid the problems of sensitivity and premature 
convergence. A unique optimal solution is obtained even when the objective functions 
are taken one at a time. 
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Chapter 1: INTRODUCTION 

Present day chemical engineering is associated with core competencies in four major 
areas: reaction engineering, transport phenomena, separations science, and computational and 
systems science. Several paradigm shifts have taken place in this discipline over the years. These 
include the introduction of mathematical modeling in its various forms (including process 
control, systems approach, etc.), the shift from unit operations to transport phenomena, the recent 
transition towards biosystems, etc. However, one facet of chemical engineers remains 
unchanged, namely that they have a responsibility of integrating the chemical engineering core 
constituencies with economic parameters so as to achieve commercial success. In this context, 
optimization of chemical processes plays a key role in chemical engineering. Optimization 
techniques have been applied to problems of industrial importance ever since the late 1 940s, and 
several excellent texts (Beveridge and Schechter, 1970; Bryson and Ho, 1969; Edgar and 
Himmelblau, 1988; Lapidus and Luus, 1967; Ray and Szekely, 1973) describe the various 
methods with examples. The complexity of the problems studied has increased as fester and 
more powerful computational hardware and software have become available. The solutions of 
even more sophisticated optimization problems will become available as complex biosystems are 
studied and get exploited commercially in the future. In this decade, particularly, several 
industrial systems have been optimized, which involve multiple objective functions and 
constraints, using a variety of mathematical techniques and robust computational algorithms. In a 
few cases, the optimal solutions have been implemented in industry with some success. 

Optimization of chemical processes has been a fascinating field of study for several 
decades. Until about 1980, virtually all problems in chemical engineering were optimized using 
single objective functions. Often, the objective function (also called the cost function) accounted 
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for the economic efficiency only, which is a scalar quantity. The reason for solving such 
relatively simple optimization problems was possibly, the limitations posed by the technology of 
computing systems. Most real world chemical engineering problems require the simultaneous 
optimization of several objectives (multiobjective optimization), which cannot be compared 
easily with each other (are non-commensurate), and so cannot be combined into a single, 
meaningful scalar objective function. Examples include reliability, safety, hazard analysis, 
control performance, environmental quality, etc., apart from the major goal of achieving 
economic efficiency. Until a few years ago, these several objective functions were combined into 
a single scalar objective function, using arbitrary weight factors, so that the problem could 
become computationally tractable. This ‘scalarization’ of what is really a vector objective 
function suffers from several drawbacks. One is that the results are sensitive to the values of the 
weighting factors used, which are difficult to assign on an a-priori basis. What is even more 
important is the less-recognized fact that there is a risk of losing some optimal solutions 
(Chankong and Haimes, 1983; Haimes, 1977). Therefore, it may be necessary that a vector form 
of all the objective functions be used in formulating and solving real-life optimization problems. 

The concept of multiobjective optimization is attributed to the economist, Pareto (1896). 
After several decades, this concept was recognized in operations research and has recently 
become popular in engineering. There have been several applications of multiobjective 
optimization in chemical engineering. Polymerization reactor optimal design and control 
problems are inevitably multiobjective in nature. Unique problems in polymerization reactors are 
encountered in the area of optimal design and control. This is because of the fact that the quality 
of the polymer is defined in abstract terms such as strength, processability, stiffness, crimp, etc. 
Though these terms can be correlated to more quantifiable properties like molecular weight 
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distribution and the composition, the definition of meaningful optimization problems always 
ends up in the formulation of a multiobjective optimization problem. 

The aim of this study is to apply multiobjective optimization on the industrial 
polyethylene terephthalate reactor. Polyethylene terephthalate (PET), commonly referred to as 
polyester, is one of the most commonly used polymers with diverse applications like synthetic 
fibres, films , bottles, engineering applications, etc. Commercially, PET is manufactured in three 
stages using continuous reactors. The first stage (esterification stage) is carried out at 
atmospheric pressure and at 270-280°C. The raw materials commonly used are a molar excess of 
ethylene glycol (EG) and either purified terephthalic acid (PTA) or dimethyl terephthalate 
(DMT). PTA and DMT have their own advantages and disadvantages. Nowadays, DMT is not 
usually preferred because of the problems in handling methanol, a by-product in the first stage. 
The use of PTA posed a problem earlier because its solubility in EG is low. But it has good 
solubility in BHET (bis-hydroxy ethylene terephthalate), a product of the first stage reactor, and 
the use of a recycle stream (or back mixing) in the first stage can overcome this drawback. PTA 
and EG are now usually processed in a series of CSTRs or a plug flow reactor with a recycle in 
the first (or esterification) stage. A polycondensation catalyst, antimony trioxide, is injected in 
small concentrations (0.03 - 0.05 percent by weight) into the oligomer stream exiting from this 
reactor. The second stage (pre-polymerization) is carried out either in one or two agitated vessels 
under reduced pressures, at about 15-30 mm Hg and 270-280°C. The degree of polymerization, 
DP, reaches a value of about 30-40 in this stage. The prepolymer so produced undergoes final 
polycondensation in a finishing (or third stage) reactor in which the pressure is maintained quite 
low at 1-2 mm Hg, at a temperature of about 280-290°C. Since the reaction mass is very viscous 
under these conditions, the finishing reactor has a special construction to enhance mass transfer 
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and the removal of the by-product, ethylene glycol, so as to drive the reaction in the forward 
direction and to give a product having a high value of DP. The finisher is usually a jacketed 
cylindrical vessel with a horizontal agitator, with large screens mounted on the latter. The 
reaction mass in the third-stage reactor is usually heated by condensing Dowtherm™ vapors in a 
jacket. The desired DP of the product (textile grade) from this reactor is above about 80. In the 
finishing reactor, a complex series of reactions accompanies the main polycondensation reaction. 
Apart from playing an important role in controlling the product quality (such as fiber color, 
dyeability and stability), these side reactions also modify the progress of polymerization. There 
is an absolute need for meeting the stringent industrial product specifications through more 
effective operation these days. For example, a small change in di-ethylene glycol end group 
concentration, [Edeg], in the polymer affects polymer crystallinity, melting point and dyeability 
of the polyester. Such polymer properties are very difficult or impossible to measure on-line. 
Thus, it is very difficult to control the polymer product quality due to delayed off-line laboratory 
measurements. Therefore, any insight into the optimal operation of the finishing stage leading to 
improved productivity and product quality control is extremely useful. 

Having understood the need for optimization of finishing reactors, one must decide on 
suitable methodologies to solve the problem. Numerous optimization techniques are available for 
solving multiobjective optimization problems (Chankong and Haimes, 1983; Srinivas and Deb, 
1995). These days, with the advent of high speed computing machines and the availability of 
powerful computational software, there has been a shift towards evolutionary techniques that 
yield comparatively much better results for complex problems with multiple solutions. One such 
technique that has rapidly gained acceptance in the recent past is Genetic Algorithm (GA). GA is 
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noted for its robustness and ability to find the global optimum where most other methods end up 
in local optima. The reasons favoring the use of GA are: 

> Objective functions can be multimodal or discontinuous. 

> Only information on the value of the objective function is required; gradient evaluation is 
not required. 

> A starting point is not needed, and the search is from a randomly generated population of 
points. 

> It is better suited to handle large problems involving many variables. 

Moreover, GA can be adapted to effectively handle real-life problems involving the optimization 
of two or more objective functions simultaneously, as applied in this study. Non-dominated 
sorting genetic algorithm (NSGA - I), a popular adaptation of GA (Srinivas and Deb, 1995), has 
been very efficient in solving multiobjective optimization problems. Indeed several problems in 
chemical engineering have been solved using NSGA - I (Garg and Gupta, 1999; Mitra et al., 
1998). 

The model developed and tuned on some industrial data [Bhaskar et al. (2000a)] was 
re-tuned in this study since some minor errors were found in the previous study. This is 
described in Chapter 2. Chapter 3 describes NSGA - I as well as a more recent adaptation 
referred to as the elitist Non-dominated Sorting Genetic Algorithm, NSGA - II [Deb and 
co workers (2000)]. Chapter 4 compares the results of the multiobjective optimization problem 
using both NSGA - 1 and NSGA - II. Finally, appropriate conclusions are made on the study and 
suggestions for future work are indicated in Chapter 5. 



Chapter 2: MODELING AND SIMULATION 


2.1 Introduction 

Poly (ethylene terephthalate) (PET) is one of the most widely used polymers in the fiber and 

engineering plastics industry. In the past two decades, PET has gained sufficient commercial 

importance to inspire a considerable amount of research in the area of modeling and simulation 

of the several reactors used in its manufacture. The early 1980s saw work in the open literature in 

this area primarily from two groups in India, ours (Kumar et al., 1982 a, b; 1983) and that of 

Ravindranath and Mashelkar (1982 a, b). After a short period of reduced interest, this area again 

started attracting considerable interest from groups led by Choi (Besnoin and Choi, 1989; Lei 

and Choi, 1990; Laubriet et al., 1991; Cheong and Choi, 1995, 1996), Ray (Hipp and Ray, 1996; 

Jacobsen and Ray, 1992a) and Doherty and Malone (Steppan et al., 1990). Several reviews have 

also appeared over these last two decades in this area (Gupta, 1983; Ravindranath and 

Mashelkar, 1986 a, b; Jacobsen and Ray, 1992b). The finishing reactor for PET poses interesting 

problems in modeling due to the coupling of the mass-transfer and kinetic effects. In the past few 

years, several attempts have been made to develop models for this reactor. These vary in their 

levels of sophistication, from the earliest and simplest models (Secor, 1969; Hovenkamp, 1971; 

Amon and Denson, 1980; Ravindranath and Mashelkar, 1982c, 1984; Ghosh et al., 1983; Gupta 

et al., 1984; Kumar et al., 1984a) to more recent ones (Laubriet et al., 1991; Cheong and Choi, 

1995, 1996; Jacobsen and Ray, 1992a; Hipp and Ray, 1996) incorporating not only the main 

reactions but also several side reactions taking place in the reactors. Secor (1969) was the first to 

use penetration theory to explain mass transfer effects in a finishing reactor. The calculations use 

a value of 1.0 for the equilibrium constant for the polycondensation reaction. Hovenkamp (1971) 

made use of the Flory-Huggins (Flory, 1953) model instead of Raoult’s law, to predict the 

6 



activity coefficients of the different vaporizing species, EG, as well as other low molecular 
wei ght compounds produced by the side reactions, namely, water (W) and diethylene glycol 
(DEG). Amon and Denson (1980) proposed a simplified analysis for the performance of wiped 
film reactors. Gupta et al. (1983) and Kumar et al. (1984a) made modifications in this model and 
extended it for predicting the poly-dispersity index of the polymer using moment-closure 
techniques (Gupta and Kumar, 1987). Ravindranath and Mashelkar (1982a) developed an 
effective flash procedure for determining the interfacial concentration of the volatile species 
along the length of the finisher. Disc-ring reactor models and axial dispersion plug flow models 
were proposed by Dietze and Kunhe (1969) and Ravindranath and Mashelkar (1982c) for 
describing the performance of these reactors. Laubriet et al. (1991) proposed a two-phase model 
for the finishing reactor for the polycondensation of PET. This model has the advantage of being 
independent of the specific reactor geometry, since it uses a single variable, ki,a, to describe the 
mass transfer effects in the reactor. Cheong and Choi (1995, 1996) proposed a multi- 
compartment model for the continuous rotating disc polycondensation reactors. This model does 
not consider the side reactions that occur in the reactor. A detailed kinetic modeling study has 
been done by Jacobsen and Ray (1992a) who solve the equations using their CAD package, 
POLYRED. Hipp and Ray (1996) have developed a dynamic model for polycondensation in 
tubular reactors. 

Bhaskar et al.(2000a) used the two phase model of Laubriet et al.(1991) and tuned the 
model-parameters so as to describe a third stage wiped film industrial PET reactor. A few 
errors were detected in the computer code used. In this study, we re-tuned ten parameters of 
Bhaskar et al. (2000a) using genetic algorithm (GA). This re-tuned model is then used to obtain 
optimal operating conditions for this reactor, using multiple objective functions and end point 
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constraints. The technique used in this study is quite general, and parallels similar work on the 
modeling of an industrial nylon-6 semibatch reactor (Wajge et al., 1994). 

2.2 Formulation 

The kinetic scheme used in the present work is presented in Table 2.1. Apart from the 
main polycondensation reaction, the side reactions considered are the formation of acetaldehyde 
(A), DEG, water (W) and the degradation of di-ester groups (Z). The reactions taking place in 
the finishing reactor are quite complex. Unfortunately, kinetic information is not available for 
many of the side reactions, and even for the main reactions is, at best, unreliable. The choice of 
the side reactions is based on the objectives of modeling and on the level of sophistication we 
wish to have in predicting the product properties. We have considered all the main reactions 
(Equations 1-9, Table 2.1). These have also been used for modeling by Bhaskar et al.(2000a). 
The concentrations of several of the other by-products are extremely low under conditions of 
interest and therefore, in all engineering applications and for computational ease, one has usually 
ignored them. For example, cyclic oligomers have been ignored by Ravindranath and Mashelkar, 
(1986a) and Laubriet et al. (1991). 

The steady-state model equations are quite similar to those presented by Laubriet et al. 
(1991) and the complete set is summarized in Table 2.2. It is assumed that the flow of the liquid 
down the reactor is plug-flow and that the vapor phase is well mixed. Also, it is assumed that 
there is no mass transfer resistance in the vapor phase. Flory-Huggins’ theory (Flory, 1953) is 
used to describe the activity coefficients of three of the four volatile species, EG, W and DEG. 
The other species, acetaldehyde, is assumed to vaporize as soon as it is produced in the liquid 
phase. Most of the details of the model development are available in the literature [Bhaskar et al. 
(2000a)] and are not repeated here. Bhaskar et al. (2000a) tuned nine parameters in the model. 
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TABLE 2.1: Kinetic scheme and kinetic parameters 


1 . Ester interchange reaction ( Main Poly condensation) 



2. Acetaldehyde Formation 



3. Diethylene Glycol Formation 




5. Degradation of Diester group 
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Table 2.1 Continued 

Kinetic parameters (Bhaskar et al. (2000))([Sb203] = 0.04 wt %) + 


Reaction 

No., i 

Type of 

reaction 

Rate constant 

ki = k i0 exp(-Ei/RT) 

Equilibrium 

constant 

k io (m 3 mol' 1 min' 1 ) 

E, (J mol' 1 ) 

1 

Reversible 

1.09xl0 3 

77441 

5.1400 

2 

Irreversible 

2.7480x1 0 7 * 

124742 


3 

Irreversible 

2.245 xlO 3 

77441 


4 

.. 

Irreversible 

8.32x1 0 4 

124742 


5 

Reversible 

1.09xl0 3 

77441 

2.871 xlO' 2 

6 

Irreversible 

8.32x1 0 4 

124742 


7 

Reversible 

2.08x1 0 3 

73673 

2.5 

8 

Reversible 

2.08x1 0 3 

73673 

j 

10.340 

9 

Irreversible 

0.6395xl0 9 * 

158230 



+ Values of Ki, K5, Kg, ko 0 , k 9o and k3 0 are those obtained in the present study, and differ from those 
of Bhaskar et al. (2000a) . 

*min' 1 
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TABLE 2.2: Complete set of model equations used in this study 


Balance equations for the liquid phase [Laubriet et aL (1991)1 


[ — 2R, — R~ — Rs* — Ra — Rr — 2 R,- 4 - R~ — RoJ 


[R 2 ^4 ■*" Rq Rfi + RgJ 


1 d[E g ] _ 

0 dz 
1 d[Eq] _ 

6 dz 

1 d[E v ] 

T~dT~ L 
1 d _EDEGL 

0 dz 

= R \ - r a - R 1 - k i a([EG] ~ [EG * ]) 

1 d[W] 


[ + Rq] 

[~r }+ r 6 ] 


6 dz 
1 d[DEG] 

0 


R 1 +R 2 -k l a([W]-[W*J) 

= R a + R, - k,a([DEGJ - [DEG*]) 
dz 4 3 / 


where : 

Ri = ki [E g ] 2 -4k 1 ' [Z] [EG] 
R 2 =k 2 [E g ] 

R 3 =k 3 [E v ][E g ] 

R4=2k4[Eg][EG] 

R 5 =k 5 [E g ] [E D EG]-4k 5 ' [Z] [DEG] 
R6=k6[Eg] 2 

R 7 =2k 7 [E a ][EG]-k 7 '[E g ][W] 

R s =k 8 [E g ][E a ]-2k 8 '[Z][W] 

R 9 =k 9 [Z] 
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Vapor-liquid equilibrium correlations: 


* 

C j = 


C 




c 


poly 




J 


x* ; C* = [EG*], [W*J, [DEG*] 


poly 

* Py . J 

x . = — 


V j J 

[E g ] + [E a ] + [E v ] + [E deg ] 


; j = EG, W, DEG 


J P°y . 

J J 

In P° = 49.703 -f 8576.7/7) -4.042 In T 

, p0 4047.606 

In P„, =18.568 
W 


In F± = 17.0326 - 
DEG 


T -33.3 
4122.52 


y,- = 


T- -122.5 

1 * 
[k^C -C )dz 

0 ■> J 


j 1 


1 


j = EG, h 


I \kf(C / ~ C i> dz + J<W + *3 fEJ[E ])dz 
j 0 7 J 0 ~ s 






1 

Y ; = exp 

J m . 

J 

V p . 

P J 
m . = — — — 

J M . 

J 


1 


1 +2i 

m . 1 

V j J 


U 


M .C . 

z — J —^~ 
J P J 
\-u 


; j = EG, W, DEG 

; j = EG, W, DEG 

; j = EG, W, DEG 


V = — 
P C 


poly 


h =a o +b o 


r \ 

U 

1 * 


u 


k 7 a = k 7 a - 
/ l ref 


o 

\ v ) 

f \a 
N 


N , 
v re / y 


DP = 


C^] + [ £ g] + [ £ v] + [ £ D £ Gl + 2 ^W 


' DEG 
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We have re-tuned all these nine parameters. Upon re-tuning, however, we found that the value 
of [E v ] 0 ut, the concentration of vinyl end groups in the product, does not match the industrial 
value. As a result, a tenth parameter, k 30 , was also tuned in the present study. The model 
equations (Table 2.2) for the reactor are a set of ordinary differential equations (ODEs). These 
are to be integrated using appropriate values of the feed conditions. These are given in Table 2.3 
[Bhaskar et al (2000a)]. The NAG subroutine, D02EJF, was used to solve these equations. The 
subroutine uses Gear’s algorithm (Gupta, 1995) with a tolerance of 10' 12 . An iterative solution is 
necessary because the liquid and vapor phases have different kinds of flow (plug and 
well-mixed). Values for the average mole fractions of the three of the four volatile species in the 
vapor phase, yj, are assumed (see Table 2.4 for the starting guess values) and then the ODE’s are 
integrated from the dimensionless value of z = 0 (feed end) to z = 1 (product end). The vapor 
phase values of yj are then re-evaluated and used as the guess values for the next iteration. This 
method of successive substitutions (Picard iteration) (Gupta, 1995) is continued till the sum of 
the squares of the yj does not change between two successive iterations by more than 1x10' . 
The program was run on a P3 (586 MHz) computer, on a Linux platform. The model provides 
values of DP and the concentrations of the hydroxyl end groups (E g ), acid end groups (E a ), di- 
ester end groups (Z), vinyl end groups (E v ), DEG end groups (Edeg), ethylene glycol (EG), water 
(W) and di-ethylene glycol (DEG), as a function of the axial position in the reactor after 
convergence is achieved. 

Four sets of industrial data were available. These are the same as used by 
Bhaskar et al. (2000a) and are given in Table 2.5. In order to obtain best-fit (tuned) values of the 
parameters in the model, an error function, I, is defined. This comprised of a sum of square 
errors between the model predicted and industrial values: 
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Table 2.3 Feed conditions [Bhaskar et al. (2000a)] 


Feed concentrations, 

Kmol/m 3 

Densities of pure liquid 
components, (kg/m 3 ) 

[E g ]f =4.0x1 O' 1 

Peg — 1 108 

[E a ] f = 2.57xl0" 3 

pw = 1000 

[Z]f =11.2 

Pdeg “1118 

[E v ]f = 1.17xl0' 3 


[EoEG]f = 0.17 


[EG]f =6.5x1 O' 3 


[W] f =4.6x1 O' 4 


[DEG] f = 4.0x1 O' 4 


DP f = 40 * 



* no unit 
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Table 2.4 Values of parameters used in this study 


Initial guess values used 
f° r yj : 

Yeg = 0.93 

yw =0.0269 

Ydeg = 0.0246 


Weighting factors used in 
eqs. 2.3, 2.4 and 2.6 

Wll = W13 = W14 = lxlO 5 

W21 = W31 =W41 = lxlO 3 


GA parameters used in this study: 


Pc = 0.9 
P m = 0.05 
Sr = 0.149 
nparam= 10 
N p = 100 

Ngen= 100 
lchrom — 320 
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( 2 . 1 ) 


7 (“) = Z W » 

Uj 



Sij,ind J 


where S u is the value of the i th property at the output (or product) end for set no. j, the subscripts, 
m and ind, represent values predicted by the model (for any set of parameters, u) and the 
industrial values, respectively, and Wij are the weight factors indicating the emphasis to be 
accorded any error of the i th property of set no. j. Genetic algorithm was used to obtain the values 
of the model parameters which minimize the error, I. The NAG subroutine D02EJF with a tolerance 
of lx 10‘ 12 was used. Several sets of parameters were tried in this study till a reasonable tuning of the 
model was accomplished. These are described in the next section. 

2.3 Results and Discussion 


The corrected computer code was tested for the reference case (see Table 2.5), using the 
values of the parameters of Bhaskar et al. (2000a). It was found that the model values of the 
product DP was higher, while for the DEG and the acid end group concentrations at the product 
end, it was lower. So, we re-tuned the nine parameters (ki, a re f, ao, bo, Kj, K5, Kg, kjo, k 9o and a) 
in the model. The objective function used was 


I = w, 


1 - 


82.00 


\ 2 


+ w. 


21 


1- 


\ J ^ DEg \ 1.) 

0.17 


+ w. 


31 


1 


[E a l, 


1.038x 10~ 3 


v r DP y 

-j 5,m 


+ w , 


1.3 


V 


82.70 


+ w, 


14 


DP. 

j 4 ,m 

v 82.30 


( 2 . 2 ) 

The parameters, kLa re f, ao, bo, Ki, K 5 , Kg, k. 20 , k 9o and a were tuned for sets 1, 3 and 4 by 
minimizing I. It was found that value of [E v ] ou t, the output concentration of vinyl end group did 
not match with the industrial value. Several simulations were performed to study the effect of the 
several rate and equilibrium constants used in this work, individually. The frequency factor (k3 0 ) 
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Table 2.5 : Industrial data and predicted values at the outlet end for the 4 data sets* 


Set# 

Operating 

Conditions 

Product Property 

Industrial 

Value 

Model predicted 

value 

1 

Reference (ref) # 

DPout 

82.00 

82.07 



[Edeg] out, kmol/m 3 

0.17 

0.1686 



[E a ] out, kmol/m 3 

1.038x 10' 3 

1.039x 10‘ 3 



[E v ] out, kmol/m 3 

7.81 17x Iff 4 

7.8252x Iff 4 

2 

Tref + 1 K 

DPout 

82.60 

82.62 

3 

P re f- 0.5 mmHg 

DPout 

82.70 

82.70 

4 

N re f+ 0.1 rpm 

DP ou t 

82.30 

82.30 


* Tuned values of Table 2.1 used to generate model-results. 

# T re f = 566.15 K; P re f = 2.0 mm Hg; [Sb 203 ] re f = 0.04 wt. %; 0* = 1.0; N* = 1.0; For the other 
sets, only one of these values is different, as indicated. 




of the rate constant, k 3 was found to affect [E v ] ou t most significantly. Hence, k 3 o was also 
included as the tenth tuning parameter. The objective function used finally is, thus, given 
by 


I = viy 


DR 


1 ,OT 


82.00 


+ w. 


21 


+ W v . 


1 -- 


DR 


3 ,m 


82.70 


J 

\ 2 


+ w, 


_ 3 l,i 

V 0-17 , 

DP. m V 


+ w 


31 


1- 


[E a ] 


1.038x10“ 


+ w. 


41 


1 -- 


[E v ] 


7.8117x10“ 


82.30 


(2.3) 


The parameters to be tuned are kLa re & ao, b 0 , Ki, K 5 , IQ, k 2 0 , k 9 0 , a and k 3 o. Table 2.6 lists the 
upper and lower bounds as well as the final tuned values (using sets 1, 3 and 4) of these 
parameters. The values of the six kinetic/equilibrium constants are also included in Table 2.6 
along with the other (original, Bhaskar et al. (2000)) values. A similar tuning of rate constants as 
well as a few other parameters have been reported by Zabisky et al. (1992), for an industrial 
LDPE reactor. The weight factors were chosen by trial and error, to achieve the best curve-fit. 
The values of the weight factors used are also given in Table 2.4. The FORTRAN 77 program 
was run on a P3 (586 MHz) computer. 

The parameters so obtained could explain the three sets of industrial results (sets 1, 3 and 4) 
used for tuning, quite well. These are then used without change , to predict the effect of a change 
in temperature (set 2). Very good agreement of the value of DP of the product for this case was 
obtained (see Table 2.5). The fact that the parameters were obtained using three sets of industrial 
data, and then the model predicted another set of industrial data well, without re-tuning of these 
parameters, indicates that the model accounts for all the physico-chemical phenomena actually 
present in the industrial reactor quite well. 
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Table 2.6 Bounds and final tuned values used in this study 
Bounds and final tuned values for the parameters describing the reactor/reactions: 


Bounds 

Tuned Values 

Tuned Values of Bhaskar et al 
(2000) 

3.4 < kL&ref ^ 3.8 

3.7670 

2.6875 

0.59 < a<, < 0.79 

0.6189 

1.0378 

2.0 < b 0 <2.1 

2.0210 

2.1838 

4.9 < Ki < 5.4 

5.1400 

6.1835 

0.027 <K 5 < 0.031 

2.871xl0’ 2 

5.143xl0' 2 

10.0 < Kg < 15.0 

10.340 

11.87 

2.4xl0 7 <k 2o < 2.8x1 0 7 

2.748x1 0 7 

4.7674x1 0 7 

5.0x l0 s < k 9o < 7.0x1 0 8 

6.395x10 s 

0.2215x10 s 

2.0 <a <4.0 

2.627 

2.6647 

l.OxlO 6 < k 3o < 4.0x1 0 6 

2.245x1 0 6 

Not Tuned (1.09xl0 6 ) 


) 
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A sensitivity analysis is now performed by varying one parameter at a time with other 
parameters fixed at their reference values. The results are shown in Table 2.7. Such a study gives 
an intuitive feel of the reactor behavior, information that is useful in studies involving reactor 
optimization and control. The effect of temperature is considered first. The concentration profiles 
of the side products and the variation of DP along the axial position, z, of the finishing reactor 
for three different temperatures (564 K, 566 K and 569 K) are presented in Figs. 2. 1-2.8. Fig. 2.1 
shows the effect of the reactor temperature on the DP of the polymer. It is observed that as the 
temperature goes up from 564 to 569 K, the value of the DP of the product increases from 80.95 
to 83.71. This is the combined effect of temperature on the rate constants, as well as on the 
equilibrium interfacial concentration, C,* (through Pj°), of the volatile species. It is assumed that 
Xi as well as ki.a re f are not much influenced by such a change in temperature. The DP of the 
polymer increases primarily because of the faster rate of the main (No. 1) reaction in Table 2.1 at 
higher temperatures, which leads, simultaneously, to a faster rate of vaporization of EG, the 
concentration of which is lower in the melt (Fig. 2.2). The degradation reaction (No. 9, Table 
2.1) is found to be strongly influenced by temperature. Both the acid end group concentration, 
[E a ], as well as the vinyl end group concentration, [E v ], increase with an increase in the 
temperature (see Figs. 2.3 and 2.4). These effects agree well with observations made in industry. 
The vinyl end groups lead to the formation of anhydrides and network structures, and so need to 
be kept low. Figs. 2. 5-2.7 show that [E g ] and the concentration of [DEG], decrease with 
increasing temperature, while concentration of [W] increases with increasing temperature. [Eg] 
decreases because the rate of the main reaction (No. 1) increases. The concentration of [DEG] is 
lower at higher temperatures because of higher rates of vaporization, as was the case with EG. 
Lower operating temperatures are preferred in industrial operations to avoid high concentrations 
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Table 2.7 Sensitivity analysis 


Concentration of various species at the reactor outlet in kniol/m 


I 


neter Value [E g ] ou t 

564 0.1135 

566 0.1095 

569 0.1038 

0.50 0.1039 

1.00 0.1057 

1.50 0.1076 

2.00 0.1095 

0.92 0.1157 

0.96 0.1125 

1.00 0.1095 

1.04 0.1067 

3.00 0.1280 

3.50 0.1177 

4.00 0.1095 

4.25 0.1060 


0.90 0.1095 

1.00 0.1095 

1.05 0.1095 

1.25 0.1094 

1.00 0.1095 

1.05 0.1088 

1.15 0.1077 

1.25 0.1067 

3.000 0.1126 

3.500 0.1104 

3.767 0.1095 


10 3 [ E a ] out 

lOWout 

0.9890 

7.3987 

1.0394 

7.8252 

1.1278 

8.5848 

1.0230 

7.9791 

1.0282 

7.9277 

1.0336 

7.8764 

1.0394 

7.8252 

1.0373 

7.6413 

1.0379 

7.7309 

1.0394 

7.8252 

1.0417 

7.9237 

0.9851 

7.1064 

1.0127 

7.4868 

1.0394 

7.8252 

1.0523 

7.9807 


1.0395 

7.8232 

1.0394 

7.8252 

1.0394 

7.8262 

1.0392 

7.8300 

1.0394 

7.8252 

1.0373 

7.8420 

1.0336 

7.8718 

1.0307 

7.8971 

1.0749 

7.6925 

1.0500 

7.7848 

1.0394 

7.8252 


[ EdegW 

10 5 [ W] out 

0.1685 

1.0427 

0.1686 

1.0881 

0.1687 

1.1678 

0.1675 

0.8457 

0.1680 

0.9236 

0.1683 

1.0044 

0.1686 

1.0881 

0.1686 

1.1338 

0.1686 

1.1098 

0.1686 

1.0881 

0.1686 

1.0685 

0.1685 

1.2088 

0.1685 

1.1411 

0.1686 

1.0881 

0.1686 

1.0650 


0.1686 

1.0895 

0.1686 

1.0881 

0.1686 

1.0873 

0.1685 

1.0845 

0.1686 

1.0881 

0.1685 

1.0575 

0.1683 

1.0050 

0.1682 

0.9622 

0.1689 

1.3226 

0.1687 

1.1570 

0.1686 

1.0881 


10 5 [DEG]out 

DPout 

1.1974 

80.95 

1.1557 

82.07 

1.0959 

83.71 

1.0759 

84.06 

1.1039 

83.37 

1.1303 

82.70 

1.1557 

82.07 

1.2208 

80.27 

1.1872 

81.19 

1.1557 

82.07 

1.1261 

82.91 

1.3508 

77.04 

1.2423 

79.76 

1.1557 

82.07 

1.1187 

83.09 


1.1562 

82.05 

1.1557 

82.07 

1.1554 

82.08 

1.1543 

82.11 

1.1557 

82.07 

1.1468 

82.29 

1.1311 

82.67 

1.1176 

83.01 

1.1905 

81.08 

1.1661 

81.77 

1.1557 

82.07 
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4.0 

0.1144 

1.0247 

7.6433 

0.1686 

1.1205 

1.2077 

80.66 

5.0 

0.1100 

1.0379 

7.8067 

0.1686 

1.0913 

1.1609 

81.92 

6.0 

0.1069 

1.0475 

7.9222 

0.1686 

1.0712 

1.1287 

82.82 

6.5 

0.1057 

1.0514 

7.9682 

0.1686 

1.0634 

1.1161 

83.17 


0.15 

0.1098 

1.0383 

7.8117 

0.1697 

1.0911 

0.6129 

81.63 

0.25 

0.1096 

1.0391 

7.8216 

0.1689 

1.0889 

1.0104 

81.95 

0.30 

0.1095 

1.0395 

7.8265 

0.1685 

1.0878 

1.2059 

82.11 

0.35 

0.1093 

1.0399 

7.8313 

0.1681 

1.0867 

1.3993 

82.26 

9.0 

0.1095 

1.0717 

7.8249 

0.1686 

1.0865 

1.1557 

82.06 

10.0 

0.1095 

1.0468 

7.8251 

0.1686 

1.0877 

1.1557 

82.07 

11.0 

0.1095 

1.0264 

7.8253 

0.1686 

1.0887 

1.1557 

82.07 

12.0 

0.1095 

1.0093 

7.8255 

0.1686 

1.0895 

1.1556 

82.08 

2.00 

0.1097 

0.9591 

7.8203 

0.1686 

1.0158 

1.1574 

82.04 

3.00 

0.1094 

1.0665 

7.8269 

0.1686 

1.1123 

1.1551 

82.08 

4.00 

0.1092 

1.1737 

7.8334 

0.1686 

1.2084 

1.1527 

82.11 

5.00 

0.1090 

1.2809 

7.8400 

0.1686 

1.3044 

1.1504 

82.14 

1.0 

0.1096 

1.0390 

13.0170 

0.1686 

1.0987 

1.1574 

81.87 

2.0 

0.1095 

1.0394 

8.5652 

0.1686 

1.0883 

1.1559 

82.04 

3.0 

0.1094 

1.0395 

6.0973 

0.1686 

1.0877 

1.1552 

82.13 

4.0 

0.1094 

1.0396 

4.6564 

0.1686 

1.0874 

1.1550 

82.18 

5.0 

0.1096 

0.9486 

6.6074 

0.1686 

1.0088 

1.1569 

82.10 

6.0 

0.1095 

1.0137 

7.4803 

0.1686 

1.0656 

1.1560 

82.08 

7.0 

0.1094 

1.0788 

8.3535 

0.1686 

1.1224 

1.1552 

82.06 

8.0 

0.1094 

1.1440 

9.2271 

0.1686 

1.1791 

1.1543 

82.03 
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Figure. 2.1 Effect of temperature on the DP of the polymer as a function of the 
dimensionless axial position in the reactor, z. The unmarked curve represents the variation 
of DP vs. z at the reference temperature of 566 K. 



Figure. 2.2 Effect of temperature on the variation of the ethylene glycol concentration, 

[EG], with z. The unmarked curve corresponds to 566 K. 
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Dimensionless position, z 


Figure. 2.3 Effect of temperature on the variation of the acid end group concentration, 
[E a ], with z. The unmarked curve corresponds to 566 K. 



Dimensionless position, z 

Figure. 2.4 Effect of temperature on the variation of the vinyl end group concentration, 
[E v ], with z. The unmarked curve corresponds to 566 K. 
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Figure. 2.5 Effect of temperature on the variation of the hydroxyl end group 
concentration, [E g ], with z. The unmarked curve corresponds to 566 K. 



Figure. 2.6 Effect of temperature on the variation of the water concentration, [W], with 
z. The unmarked curve corresponds to 566 K. 
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Figure. 2.7 Effect of temperature on the variation of the DEG concentration, [DEG], 
with z. The unmarked curve corresponds to 566 K. 



Figure. 2.8 Effect of temperature on the variation of the DEG end groups concentration, 
[Edeg]» with z. The unmarked curve corresponds to 566 K. 


26 





of [E a ], which is known to affect the hydrolytic stability of the polymer product. High 
concentrations of [E a ] lead to problems during the spinning and drawing stages in fiber 
manufacture due to the use of moist air and water in these operations. Lower temperatures also 
favor lower values of [E v ]. Fig. 2.8 shows that increasing temperatures leads to higher 
concentrations of Edeg at the outlet of the reactor. This is in contrast to the trend observed for the 
concentration of the volatile species, DEG, the concentration of which is far lower due to its 
vaporization. The increase of [Edeg] with increasing temperatures is, indeed, observed for the 
industrial reactor being studied, and is opposite to the effect observed by Laubriet et al. (1991). 
The concentration of Edeg in the final polymer determines its degree of crystallinity and hence 
its dyeability, and it is important to have a model, which predicts the right trends. Table 2.7 
summarizes the effects of varying the different operating conditions and parameters on the 
values of the concentrations of the various side products, as well as the DP at the outlet of the 
reactor. It is observed that the pressure in the reactor has a strong influence on the DP of the 
product. The DP of the polymer at the reactor outlet increases from about 82.07 to 84.06 as the 
pressure decreases from 2.0 to 0.5 mm Hg. Decreasing pressures facilitate the removal of EG 
from the bulk melt phase and enhance the growth of the polymer chains. The concentration of 
DEG in the product decreases as the pressure is lowered. This effect is similar to what happens 
with increasing temperatures. The value of [E a ] is not affected much by changing pressures, in 
sharp contrast to the behavior with increasing temperatures. The effect of decreasing pressure on 
[E v ] in the product is, again, much less than the effect of temperature. It is clear that lower 
pressures are used to produce higher molecular weight polymer since they compensate the 
effects of lower temperatures that are required to keep the other product properties within 
prescribed limits. These trends would be useful in optimization studies in the future. 
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Table 2.7 shows that the DP of the final product increases as the residence time is 
increased. This is associated with a faster depletion of hydroxyl end groups. The other side 
product concentrations are not affected much by a change in the residence time. Again, one 
would like to minimize the residence time in any optimization study, while producing polymer 
having the same values of DP, thereby increasing the throughput. Table 2.7 shows that while ao, 
the value of xi at the feed end of the reactor does not influence the DP of the product much, b 0 
does, indeed, have some effect. b 0 is also seen to affect the value of [Edeg]- Though the timed 
values of a<, and b 0 were able to predict the industrial results fairly well, the sensitivity of DP and 
[Edeg] to b 0 suggests the importance of obtaining better correlations for 

An increase in ki,a rc f leads to increasing values of DP, while it decreases [E a ], both being 
desirable objectives. The value of [Edeg] decreases with an increase in kLa re f, which may not be 
desirable. An increase in ki.a could be achieved by a better design of the agitator, an increase in 
the speed of the agitator, or by altering the level of the melt in the reactor. The first of these 
cannot be achieved except at the design stage. Lowering the level of the melt in the finisher 
simultaneously leads to an increase in the residence time, which counteracts the advantages 
attained because of lower rates of production. The manipulation of the agitator speed, N, is, 
therefore, a more preferred control-action in industry, particularly if the reduction of pressure is 
not sufficient to achieve an objective. In fact, a modification of pressure and/or the speed of the 
agitator are most often used for negating the effects of unplanned disturbances affecting the DP 
of the product in actual reactor operations. 

An increase in the catalyst concentration from 0.03 or 0.0425% leads to a significant 
increase in both the DP of the product and vinyl end group concentration, the first desirable and 
second one undesirable attribute. Indeed, changes in the catalyst concentration and temperature 
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are used as the final recourse in industrial practice when changes in both pressure and the speed 
of the agitator are unable to negate the effect of disturbances on the DP of the polymer product. 
Unfortunately, the time lags associated with both these corrective control actions (catalyst 
concentration and temperature) are large, and this is why these actions are used only when the 
other means (pressure and the speed of the agitator) have failed. 

The catalyst concentration also affects the color of the polymer produced. Barkley and 
Predmore (1970) and Chimaru et al. (1973) have suggested that the Sb J+ in the SbzCb catalyst 
gets reduced to metallic Sb° during the reaction. David and Lawrence (1995) have found 
experimentally that there is a non-linear dependence of the L-color (a measure of the luminosity 
or the whiteness) of the polymer on the concentration of Sb°. In fact, increasing SbiOs 
concentrations leads to an undesirable greying of the polymer (lower L-color). Since our kinetic 
model does not incorporate this phenomenon, we should constrain the concentration of the 
catalyst so as to achieve a balance between the L-color and the two other properties of the 
polymer (DP and [E v ]), which are strongly influenced, by the catalyst concentration. 

Since several rate and equilibrium constants have been ‘re-tuned’ in this study to achieve a 

good fit, it is important to see the sensitivity of the results to these parameters. Table 2.7 shows 

that of the four equilibrium constants present in the reaction scheme, the DP of the polymer 

product is more significantly influenced by Ki and K 5 . An increase in K 5 also leads to a 

significant decrease in [Edeg]- Indeed, K 5 was used as a tuning parameter in this study, primarily 

to achieve a better curve-fit of the industrial [Edeg] value. An increase in leads to a desirable 

lowering of the acid end group concentration. Table 2.7 also shows that increasing k 2o leads to an 

increase of the acid end group concentration, and an increase in k 3 0 leads to a significant 

decrease in [E v ] and a small increase in DP of the polymer, while an increase in k 9o results in a 

decrease of DP but a substantial increase in [E a ], [E v ] and [W], which is highly undesirable. It is 
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clear from this sensitivity study that one needs to get better estimates for the rate and equilibrium 
constants using experimental data under more ideal reaction conditions, rather than from 
industrial data where several additional physical phenomena are present. It may be added that 
tuning of the rate and equilibrium constants to explain industrial data has been used for other 
polymerization systems as well, e.g., LDPE (Zabisky et al., 1992) and nylon 6 (Wajge et al., 
1994). 

2.4 Summary 

An improved and re-tuned model for the finishing reactor in PET manufacture is 
developed. Optimal values of the parameters of the model are estimated using industrial data 
available with us. Data involving commercially important product properties like DP and the 
concentrations of the DEG and acid end groups have been used. The parameters tuned are kLa re f, 
the Flory-Huggins’ interaction parameter, xi, the equilibrium constants, Ki, K 5 and Kg, the 
frequency factors of three-rate constants, ko, k 3 and k 9> and a, a parameter associated with the 
effect of the rpm of the agitator on kLa. These parameters are tuned using three sets of industrial 
data. The tuned parameters are then used unchanged to predict the concentration of the various 
side products and the DP for a different operating temperature. It is found that the model predicts 
these values very well. A sensitivity analysis is then performed using the tuned values of the 
model parameters, and the effect of several operating variables as well as parameters on the 
important physical properties is studied. Conflicting trends are observed from the model and 
insights developed can now be used to formulate and obtain meaningful solutions to problems on 
the optimization as well as control of this industrial reactor. This study brings out the importance 
of having accurate kinetic data, as well as of values of xi and kLa. 
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Chapter 3: INTRODUCTION TO NSGA - 1 AND NSGA -II 


3.1 Genetic Algorithm (GA) 

GA is a search technique developed by Holland (1975) that mimics the process of natural 
selection and natural genetics. In this algorithm, we code a set of values of the decision variables 
(a solution, x) in terms of a ‘string (or chromosome)’ of binary numbers, generated using random 
numbers. A ‘population (gene pool)’ of such binary strings is first generated. Each chromosome 
is then mapped into a set of real values of the decision variables, using the upper and lower 
bounds of each of these. This ensures that the decision variables lie within their bounds. Then, a 
model of the process is used to provide values of the objective function for each chromosome. 
The value of the objective function of any chromosome reflects its ‘fitness’. The Darwinian 
principle of ‘ survival of the fittest’ is used to generate a new and improved gene pool (new 
generation). This is done by preparing a ‘mating pool’, comprising of copies of chromosomes, 
the number of copies of any chromosome being proportional to its fitness (Darwin’s principle). 
Pairs of chromosomes are then selected randomly, and pairs of daughter chromosomes generated 
using operations similar to those in genetic reproduction. The gene pool evolves, with the fitness 
improving over the generations. 

GA is noted for its robustness. This algorithm is superior to traditional optimization 
algorithms in many aspects, and has become quite popular in recent years. It is better than 
calculus-based methods (both direct and indirect methods) that generally seek out the local 
optimum, and which may miss the global optimum. Most of the older techniques require values 
of the derivatives of the objective functions, and quite often, numerical approximations of the 
derivatives are used for optimization. In most real-life problems, the existence of derivatives is 
questionable and often, the functions are discontinuous, multi-modal and noisy. In such cases, 
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calculus-based methods fail. Enumerative schemes, which are based on the point-by-point 
comparison of the values of the objective function in a discretized infinite (or even a finite) 
search space, are inefficient for large problems since the search space is often, too large. Random 
search techniques, too, suffer from a similar disadvantage since they work like enumerative 
techniques in the long run. GA is superior to these techniques since it is conceptually different 
from these traditional algorithms in several respects. It uses a population of several points 
simultaneously, as well as works with probabilistic (instead of deterministic) operators. In 
addition, GA uses information on the objective function and not its derivatives, nor does it 
require any other auxiliary knowledge. 

Three common operators are used in GA [called simple GA (SGA), to distinguish it from 
its various adaptations] to obtain an improved (next) generation of chromosomes. These are 
referred to as reproduction, cross-over and mutation. Reproduction, as described above, is the 
generation of the mating pool, where the chromosomes are copied probabilistically, based on 
their fitness values. Then, a pair of daughter chromosomes are produced by selecting a cross- 
over site (chosen randomly) and exchanging the two parts of the pair of parent chromosomes 
(selected randomly from the mating pool), as illustrated below for two chromosomes carrying 
information about three decision variables, each represented by four binary digits: 

1001 lilOO 0111 1001 $91 1100 

I I 

I I 

ooii oiioi iioo => ooiibiooom (3.i) 

I I 

Parent chromosomes Daughter chromosomes 

In Eq. 3.1, the crossover site for this pair of parent chromosomes is taken just after the fifth 
binary digit. It is hoped that the daughter strings are superior. If they are worse than the parent 
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chromosomes, they will slowly die a natural death over the next few generations (the Darwinian 
principle at work). 

The mutation operator is required for the following reason. In Eq. 3.1, let us assume that 
all the chromosomes in the gene pool have a 0 at the second position. There is a finite probability 
of this happening, since the generation of the binary numbers is done probabilistically. It is 
obvious that cross-over will never be able to generate chromosomes with a 1 at this position, and 
if it so happens that the optimum is, indeed, located at a point described by such a chromosome, 
GA will be unable to reach this solution. The mutation operator looks at each binary digit in 
every daughter chromosome in the gene pool, and transforms a 0 into a 1 (or vice versa) with a 
small probability. In effect, it moves the chromosome locally in the x-space, to create a better 
chromosome. The entire process is repeated till some termination criterion is met (the specified 
maximum number of generations is attained, or the improvements in the values of the objective 
functions become lower than a specified tolerance). A more elaborate description of SGA is 
available in Holland, 1975, Goldberg, 1989, and Deb, 1995. 

3.2 NSGA - I (Nondominated Sorting Genetic Algorithm - I) 

SGA has been extended to solve problems involving multiple objectives. Since GA uses a 
population of points, it seems very natural to use GA for such problems to capture a number of 
solutions simultaneously. Schaffer (1984) was the first to apply an adapted Vector Evaluated 
Genetic Algorithms (VEGA) to solve multiobjective optimization problems. Though it was 
simple to implement, the method has a bias towards some Pareto-optimal solutions (Goldberg, 
1989; Srinivas and Deb, 1995). In order to overcome this problem of bias with some of the 
optimal solutions, Goldberg proposed a non-dominated sorting procedure. The idea was 
implemented in different ways by Fonseca and Fleming (1993), Horn et al. (1994) and Srinivas 
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and Deb (1995). The algorithm implemented by Srinivas and Deb (1995) is called the non- 
dominated sorting genetic algorithm (NSGA - 1). It is to be noted that NSGA - 1 overcomes the 
pitfalls of the previous two techniques (Fonseca and Fleming, 1993, Horn et al. 1994). Fonseca 
and Fleming (1998 a,b) have recently extended their algorithm, applying the principles of niche 
and sharing. 

NSGA - I uses a ranking selection method to emphasize the good points and a niche 
method to create diversity in the population without losing a stable sub-population of good 
points. Principally, NSGA - I differs from SGA only in the way of selection. A check for non- 
dominance is first carried out among all the chromosomes in the gene pool before reproduction 
is performed. All the non-dominated chromosomes from the entire population are first identified 
and assigned a front number (Front No. = 1). These non-dominated chromosomes are then 
assigned a dummy fitness value (which is usually the number of chromosomes, N p , but could be 
any other arbitrarily selected large value instead). The dummy fitness value of any chromosome 
in this front is then modified according to a sharing procedure (Goldberg and Richardson, 1987; 
Deb, 1989; Deb and Goldberg, 1991) by dividing it by the niche count of the chromosome. The 
niche count of a chromosome represents the number of neighbors around it, with distant 
neighbors contributing less than those nearby. The niche count, thus, gives an idea of how 
crowded the chromosomes are in the x-space. This is obtained, say, for the i th chromosome, by 
computing its distance, dq, from another, say, j th chromosome in the x-space, and using a sharing 
function, Sh, as given below 


Sh(d ij ) = 1- 



\ ® share J 


If dij < ctshare ; 0, otherwise 


(3.2) 
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In Equation 3.2, cj s harc, a computational parameter, is the maximum distance allowed between 
two chromosomes to qualify as neighbors. Obviously, if dy is larger than a S hare, its contribution to 
Sh is zero (the j th chromosome is then not considered to be a neighbor of the i *), while if dy = 0, 
its contribution to Sh is 1 . For intermediate values of the distance between the two chromosomes, 
Sh lies between 0 and 1. Thus, by summing up Sh (dy) for all values of] in any front comprising 
of non-dominated chromosomes, one can get an idea of how crowded the i th chromosome really 
is. This summation is referred to as the niche count of chromosome i. The shared fitness value of 
chromosome i is the ratio of the common dummy fitness value assigned initially to each member 
of the front, and its niche count. Use of the shared fitness value for reproduction, thus, helps 
spread out the chromosomes in the front since crowded chromosomes are assigned lower fitness 
values. This procedure is repeated for all the members of the first front. Once this is done, these 
chromosomes are temporarily removed from consideration, and all the remaining ones are tested 
for non-dominance. The non-dominated chromosomes in this round are classified into the next 
front (Front No. = 2). These are all assigned a dummy fitness value that is a bit lower than the 
lowest shared fitness value of the previous front. Sharing is performed thereafter. This procedure 
is continued till all the chromosomes in the gene pool are assigned shared fitness values. The 
usual operations of reproduction, cross-over and mutation are now performed. It is clear that the 
non-dominated members of the first front that have fewer neighbors, will get the highest 
representation in the mating pool. Members of later fronts, which are dominated, will get lower 
representations (they are still assigned some low fitness values, rather than ‘killed’, in order to 
maintain the diversity of the gene pool). Sharing forces the chromosomes to be spread out in the 
x-space. The population is found to converge very rapidly to the Pareto set. It is to be noted that 
any number of objectives (both minimiz ation and maximization problems) can be solved using 
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this procedure. A flowchart describing this technique is presented in Figure 3.1 (Mitra et al., 
1998; Garg and Gupta, 1999). Additional details about the algorithm, its efficiency over other 
techniques and some comments on the choice of the computational parameters to be used in 
NSGA - 1, are described in Srinivas and Deb, 1995 and Deb, 1989. 

3.3 NSGA - II {Elitist Non-dominated Sorting Genetic Algorithms - II) 

Over the years, the main criticisms of the NSGA - 1 proposed by Srinivas and Deb(1995) 
have been : 

High computational complexity of non-dominated sorting: The currently-used non-dominated 
sorting algorithm has a computational complexity of 0(MN) (where M is the number of 
objectives and N is the population size). This makes NSGA-I a computationally intensive 
algorithm for large population sizes. This large complexity arises because of the complexity 
involved in the non-dominated sorting procedure in every generation. 

Lack of elitism: Recent results show clearly that elitism can speed up the performance of GA 
significantly. Also, it can help preventing loss of good solutions once they are found. 

Need for specifying the sharing parameter, a s hare- Traditional mechanisms of ensuring diversity 
in a population so as to get a wide variety of equivalent solutions have relied mostly on the 
concept of sharing. The main problem with sharing is that it requires the specification of a 
sharing parameter (a s h are )- Though there has been some work on dynamic sizing of the sharing 
parameter, a parameter-less diversity preservation mechanism is desirable. An improved version 
of NSGA, called NSGA - II {Elitist non-dominated sorting genetic algorithm) alleviates all the 
above three difficulties. A number of different modules that form parts of NSGA — II are 
presented below: 

A. A Fast Non-dominated Sorting Approach 
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In NSGA — I, in order to sort a population of size N according to the level of non- 
domination, each solution must be compared with every other solution in the population 
to find if it is dominated. This requires O(MN) comparisons for each solution, where M is 
the number of objectives. When this process is continued to find the members of the first 
non-dominated class for all population members, the total complexity is 0(MN 2 ). At this 
stage, all individuals in the first non-dominated front are found. In order to find 
individuals of the next front, the solutions of the first front are temporarily discounted 
and the above procedure is performed again. The procedure is repeated to find subsequent 
fronts. As can be seen, the worst case (when there exists only one solution in each front) 
the complexity of this algorithm without any book-keeping is OfMN 3 ). In the following, 
we describe a fast non-dominated sorting approach which will require at most 0(MN 2 ) 
computations. 

This approach ensures a better book-keeping strategy to make it a faster 
algorithm. In this approach, every solution from the population is checked with a 
partially filled population for domination. To start with, the first solution from the 
population is kept in a set P . Thereafter, each solution, p (the second solution onwards) is 
compared with all members of the set P , one by one. If the solution, p, dominates any 
member q of P , then solution q is removed from P . In this way, non-members of the 
non-dominated front get deleted from P. Otherwise, if solution/? is dominated by any 
member of P , solution p is ignored. If solution p is not dominated by any member of P , 
it is entered in P . This is how the set P grows with non-dominated solutions. When all 
solutions of the population have been checked, the remaining members of P constitute 
the non-dominated set. 
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Here, we observe that the second population member is compared with only one 
solution of P , the third solution with at most two solutions of P , and so on. This requires 
a maximum of 0(N 2 ) domination checks. Since each domination check requires M 
function value comparisons, the maximum complexity of this approach to find the first 
non-dominated front is also 0(MN 2 ). 

To find other fronts, the members of P will be discounted from P and the above 
procedure is repeated. 

3. Diversity Preservation 

NSGA - I uses the well-known sharing function approach, which has been found 
to maintain sustainable diversity in a population with appropriate setting of its associated 
parameters. The sharing function method involves a sharing parameter, (Tshare, which sets 
the extent of sharing desired in a problem. This parameter is related to the distance metric 
chosen to calculate the proximity measure between two population members. The 
parameter, a S harc, denotes the largest value of that distance metric within which any two 
solutions share each other’s fitness. The user usually sets this parameter. There are two 
difficulties with this sharing approach: 

1 . The performance of the sharing function method in maintaining a spread of 
solutions largely depends on the chosen value of a sh a re . 

2. Since each solution must be compared with all other solutions in the 
population, the overall complexity of the sharing function approach is 

G(N 2 ). 

In NSGA - II , the sharing function approach is replaced with a crowded comparison 

approach which eliminates both the above difficulties to some extent. This approach 
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does not require any user-defined parameter for maintaining diversity among 
population members. Also, the suggested approach is better in terms of 
computational complexity. To describe this approach, we first define a density 
estimation metric and then present the crowded comparison operator. 

B.l Density Estimation 

To get an estimate of the density of solutions surrounding a particular solution in the 
population, we calculate the average distance of two points on either side of this point along each 
of the objectives. This quantity, instance-, serves as an estimate of the size of the largest cuboid 
enclosing any point, /, without including any other point in the population (we call this the 
crowding distance). 

The crowding distance computation requires sorting of the population according to the 
value of each objective function in their ascending order of magnitude. Thereafter, for each 
objective function, the boundary solutions (solutions with smallest and largest function values) 
are assigned an infinite distance value. All other intermediate solutions are assigned a distance 
value equal to the absolute difference in the function values of two adjacent solutions. This 
calculation is continued with other objective functions. The overall crowding distance value is 
calculated as the sum of individual distance values corresponding to each objective. The 
complexity of this procedure is governed by the sorting algorithm. Since M independent sorting 
of at most N solutions (when all population members are in one front) are involved, the above 
algorithm has a computational complexity of 0(MN logN). 

After all members of the population in the set are assigned a distance metric, we can 
compare two solutions for their extent of proximity with other solutions. A solution with a 
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smaller value of this distance measure is, in some sense, more crowded than other solutions. This 
is exactly what we compare in the proposed crowded comparison operator, described below. 

B.2 Crowded Comparison Operator 

The crowded comparison operator (-< n ) guides the selection process at the various stages 
of the algorithm towards a uniformly spread-out Pareto -optimal front. Let us assume that every 
individual i in the population has two attributes: 

1. non-domination rank (i ra nk), and 

2. crowding distance (instance). 

We now define a partial order, < n , as : 

i -^n j it (irank jrank) OF [(i r ank ~ ' jrank) 3Ild (^distance - > j distance)] 

That is, between two solutions with differing non-domination ranks we prefer the 
solution with the lower (better) rank. Otherwise, if both solutions belong to the same front, then 
we prefer the solution which is located in a lesser-crowded region. 

With these three new innovations - a fast non-dominated sorting procedure, a fast 
crowded distance estimation procedure, and a simple crowded comparison operator, we are now 
ready to describe NSGA -II. 

C. The Main Loop 

Initially, a random parent population, Po, is created. The population is sorted based on 

non-domination. Each solution is assigned a fitness (or rank) equal to its non-domination 

level (1 is the best level, 2 is the next-best level, and so on). Thus, minimization of fitness 
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is assumed. At first, the usual binary tournament selection, recombination, and mutation 
operators are used to create a child population Q 0 of size N. Since elitism is introduced by 
comparing the current population with the previously found best non-dominated 
solutions, the procedure is different after the initial generation. A flowchart of a 
generation of NSGA - II is given in Fig. 3.2. 

First, a combined population, R t = P t u Q t , is formed. The population, Rt, will be of size 2N. 
Then, the population R t is sorted according to non-domination. Since all previous and current 
population members are include in R t , elitism is ensured. Now, solutions belonging to the best 
non-dominated set, / 1 , are the best solutions in the combined population and must be emphasized 
more than any other solution in the combined population. If the size of f\ is smaller than N, we 
definitely choose all members of the set, f\, for the new population, P t+ i. The remaining 
members of the population. Pm, are chosen from subsequent non-dominated fronts in the order 
of their ranking. Thus, solutions from the set fz are chosen next, followed by solutions from the 
set fs, and so on. This procedure is continued till no more sets can be accommodated. Let us say 
that the set, f\, is the last non-dominated set beyond which no other set can be accommodated. In 
general, the count of solutions in all sets from f\ to // would be larger than the population size. 
To choose exactly N population members, we sort the solutions of the last front using the 
crowded comparison operator, -< n , in the descending order and choose the best solutions needed 
to fill all the population slots. The new population, P t +i, of size N is now used for selection, 
crossover and mutation to create a new population, Qt+i, of size N. It is important to note that we 
use a binary tournament selection operator but the selection criterion is now based on the 
crowded comparison operator, -< n . Since this operator requires both the rank and crowded 
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Figure 3.2 : Flowchart of a generation of NSGA-II 


Rt(2N) - P t (N) w Q t (N) 
[Combine parent and 
child population] 


f = fast-nondominated-sort 
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i = i + 1 


Qtti = make-new-pop (Pm) 
[Use selection, crossover 
and mutation to create a 
new population Q t +i] 


t = t + 1 
increment the 
generation 
counter 















distance of each solution in the population, we calculate these quantities while forming the 
population, P t n- 

Let us now look at the complexity of one iteration of the entire algorithm. The basic 
operations and their worst-case complexities are as follows: 

1 . Non-do minated sorting is 0[M(2N) 2 ], 

2. Crowding distance assignment is [M(2N)log (2N)], and 

3. Sorting on is 0[2Nlog(2N)]. 

As can be seen, the overall complexity of the above algorithm is 0(MN 2 ), which is governed by 
the non-dominated sorting part of the algorithm. 

The diversity among non-dominated solutions is introduced by using the crowding 
comparison procedure, which is used in the tournament selection and during the population 
reduction phase. Since solutions compete with their crowding distance (a measure of density of 
solutions in the neighborhood), no extra niche parameter (such as cr S hare needed in NSGA - I) is 
required here. Although the crowding distance is calculated in the objective function space, it can 
also be implemented in the parameter space, if so desired. 

3.3 Summary 

In this section, we discussed the use of genetic algorithms for solving multiobjective 
optimization problems. Two adaptations of the genetic algorithm, NSGA - I and NSGA - II, 
were discussed and the flowchart describing each is presented. In the next section, multiobjective 
optimization of the finishing reactor under consideration is discussed. 
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Chapter 4: MULTIOBJECTIVE OPTIMIZATION OF THE 

FINISHING REACTOR 


4.1 Introduction 

Process industries aim at maximizing their production capacities while simultaneously 
maintaining the product quality. Usually, there exists a trade-off between these two 
requirements. This is particularly true in the manufacture of polymers where the properties of the 
product are crucial and reactors have to be operated under conditions which yield products that 
satisfy relatively narrow specifications. At the same time, the operating variables must be at their 
optimal conditions to maximize the throughput. It is well established that the average molecular 
weight of the polymer produced determines several important physical properties of the material, 
e.g., strength, impact resistance, etc. In addition, the concentrations of a few side products need 
to be below low limits or between narrow limits to ensure some other properties like color, 
dyeability, etc., to lie within specifications. Thus, the design and operation of polymerization 
reactors require optimization using multiple objective functions and constraints, which are often 
conflicting. A multiobjective optimization problem was formulated for the finishing reactor in 
the polyethylene terephthalate manufacture and the problem was solved using two adaptations of 
GA, NSGA - 1 and NSGA - II. 

4.2 Problem Formulation 

The two important objective functions for this system are the minimization of the acid 
and the vinyl end group concentrations, [E a ] 0 ut and [E v ] 0 ut, respectively, in the product of the 
finisher. The acid end group makes the polymer susceptible to hydrolysis [Besnoin and Choi 
(1989)] during the downstream operations and leads to breakage of the filaments during 
spinning, where the humidity is very high. The vinyl end groups have been shown to be 



responsible for the coloration of PET [Besnoin and Choi (1989), James and Packer (1995)] 
because of reactions not well understood right now, and not included in Table 2.1. Hence, the 
minimization of these two end groups improves the quality of the polymer product. The 
reduction of [E a ] 0 ut simultaneously increases the rate of polymerization of the acid end group 
catalyzed polycondensation reaction and helps maximize the throughput (this catalytic effect is 
not directly incorporated in the model, however). The important end-point constraint is to 
produce polymer having a desired value of the degree of polymerization (DP). Though the 
diethylene glycol (DEG) end groups affect the crystallinity and hence the melting point of the 
polymer unfavorably, they do improve the dyeability of the fiber. Therefore, it is preferred to 
have a certain allowable range for the concentration of DEG end groups in the fiber-grade 
polyesters. An inequality constraint is, therefore, imposed for the DEG end group concentration 
in the product, as per industrial requirements. In addition, a further inequality constraint on the 
maximum allowable limit for the acid end group concentration is imposed to ensure that it is not 
only minimized, but also lies below an upper limit. Five decision variables are used for 
optimization in this study. These are: the reactor pressure (P), temperature (isothermal, T), 
catalyst concentration ([SbaCb]), residence time of the polymeric reaction mass inside the reactor 
(0) and the speed of the wiped-film agitator (N). All of these variables can easily be changed in 
any industrial, wiped film reactor for PET manufacture, including the one being studied, and are, 
therefore, used to obtain the best optimal operating conditions. 

The multiobjective function optimization problem described above and studied in this 
work is, thus, described mathematically by [Bhaskar et al. (2000b)]: 

Min / (P, T, [SbaOs], 0*, N*) = Pi, hf = [[E a ] 0 ut, [E v ]out] T (4. la) 

subject to (s. t.) 
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DPou. = DPd 

(4.1b) 

[Ea]out< 1.038 X 1 (T kmol/W 

(4.1c) 

0. 1 660 < [EdhgW < 0.17 kmol/m 3 

(4. Id) 

cLv/dz = f(x,u); x(z=0 ) = x 0 

(4. le) 

Bounds: 


0.05 < P < 0.266 kPa (0.4 < P < 2.0 mm Hg) 

(4.1f) 

564 < T < 570 K 

(4-lg) 

0.03 < [SbiO.i] < 0.045 wt% 

(4-lh) 

0.90 <0* < 1.06 

(4-1 i) 

0.93 < N* < 1.05 

(4-lj) 


where the subscripts ‘out’ and ‘d’ refer to the values at the outlet of the reactor and the desired 
values of the product property, respectively. The variables, 0* and N*, represent dimensionless 
values, 0/0 re r and N/N re f, where 0 rc f and N re f are the values being used currently in the industrial 
reactor being studied. These two values are confidential and are not being provided here due to 
proprietary reasons. Meaningful bounds have been chosen on the five decision variables, u, 
based on industrial practice. These are given in Eqs. 4.1f-j. The mathematical model for the 
present study, the parameters and the other variables and details have already been presented in 
Chapter 2. In general, the state variable equations describing the reactor can be written in the 
form 

dx/dz =f(x,u); x(z=0)=x o (4.2) 

where x is the vector of state variables, given by 

x = [ [EJ, [EJ, [Z], [E v ], [Edeg], [EG], [W], [DEG] ] T (4.3) 
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and u is the vector of decision variables given in Eq. 4 . 1 a. 

The end-point constraint on DP in Eq. 4.1b is incorporated in the form of a penalty 
Unction (with w 2 = 10 ,s ). in the two objective functions as follows: 


Min /,* = w'i U'Uou, + w 2 


Dll 


DP, 


+ P. 


(4.4a) 


d J 


Min I* = w, V\\ na + w 2 


DP„. 


V 


DP 


+ P 


(4.4b) 


J 


An additional large ‘penalty’ value, P c (used as 10 4 in this study), is added to the objective 
functions, //* and />* if either of the two inequality constraints (Eqs. 4.1c and 4. Id), the 
maximum limit of the acid end group concentration at the outlet of the reactor, [E a ] ou t, and the 
upper and lower bounds on the DEG end group concentration at the outlet of the reactor, 
[Edi;g]oui, are violated. If these constraints are satisfied, then the value of P e in Eq. 4.4 is taken as 
zero. This ensures that bad chromosomes in the genetic algorithm used for optimization do not 
get reproduced in the successive generations even if several chromosomes violating these 
constraints exist in the initial population (i.e., chromosomes violating Eqs. 4.1c and 4.1d die out 
rapidly). Also, a very large value of P c is added if the concentrations of three of the vaporizing 
species, EG, W and DEG in the vapor phase are higher than their respective vapor-liquid 


interfacial equilibrium concentrations, EG*, W* and DEG*, at any axial location. This ensures 

vaporization (rather than condensation) of these species at all axial locations. Minimization of 

Ii* and I 2 * leads to a decrease in the acid end group concentration and the vinyl end group 

concentration in the product, while simultaneously giving preference to solutions satisfying the 

several requirements discussed above. The use of a penalty to ‘kill’ chromosomes when 

important physical constraints are violated is one of the several adaptations of the basic 

algorithm used in this study [Bhaskar et al. (2000)]. Use of wi and W 3 as 10 leads to values of 
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Wi [Ea]out &nd w^fEvjout ot the order of unity. These two objective functions are conflicting in 
nature [Bhaskar et al. (2000b)]. Use of these is, thus, expected to lead to Pareto optimal 
solutions. 

4.3 Results and Discussion 

The solution of the multiobjective optimization problem described in Eq. 4.1 is obtained 
using both NSGA — 1 and NSGA -II, using values of the computational parameters and bounds 
of the decision variables, as given in Table 4.1. The value of DP d was taken as 82.00, the value 
corresponding to the reference case for the industrial reactor simulated earlier. Such a study 
would help explore if changes in the operating variables could improve the performance of the 
industrial reactor. Figs. 4. 1 and 4.2 show all the feasible solutions (i.e., those satisfying the end 
point constraints in Eqs. 4. 1 b-d) at different values of the generation number, N g , when NSGA - 
II is used. In the initial generations, the feasible points are quite scattered in the [E a ] 0 ut - [E v ] 0 ut 
space, but as the chromosomes evolve over the generations, the feasible points move towards a 
unique point. There are 100, almost indistinguishable points in Fig. 4.3 (N g = 100). The unique 
solution of our multiobjective function optimization problem with end-point constraints (using 
NSGA - II) is given by 

[E a ]out = 9.58 x 10' 4 kmol/m 3 

[E v ]om = 7.392 x 1 0 4 kmol/m 3 

DP out = 82.000 ( 4 - 5 ) 

Table 4.2 shows the values of the five decision variables corresponding to this optimal point 

(reference run), as well as corresponding values of several other product characteristics under the 
optimal operating condition. 
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Table 4.1 : Computational variables used in this study 


w eighmUkcmm 

w, = W 3 = 1 0 4 ; W2= 10 s ; P, = 1 0 4 


QA Parameters: 

N g en =100 N P =10 ° 

p c =0.90 p m =0.001 

Wom = 1 60 bits Sr = 0.32 1 9 (NSGA - II), 0.6237 (NSGA - 1) 
Initial ( - niCSS values us ed for vy 


Yeg : 0.93 
y w : 0.0269 
Ydlx; : 0.0246 

Variables/units Lower bound Upper bound Reference value 


P, kPa (mm Hg) 

0.0532(0.4) 

0.266(2.0) 

0.266(2.0) 

T(K) 

564.00 

570.00 

564.15 

[Sb 2 0 3 ] (wt.%) 

0.03 

0.045 

0.04 

0* 

0.90 

1.06 

1.00 

N* 

0.93 

1.05 

1.00 
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Table 4.2: Optimal operating conditions and product properties * 


Operating 

conditions 

Optimal 
(NSGA - II) 

Optimal 
(NSGA -I) 

Industrial 

Bhaskar et. 
al. 

P, kPa (mm Hg) 

0.0754 

(0.5674) 

0.1787 

(1.3441) 

0.266 (2.0) 

0.185 (1.39) 

T, K 

564.05 

564.11 

566.0 

564.02 

10~[Sb 2 o7]7~ 

wt% 

3.67 

4.02 

4.0 

4.26 

8* 

1.0350 

0.9940 

1.0 

0.9672 

N* 

0.9929 

1.0201 

1.0 

1.0 

Product 

properties 





DPout 

82.00 

82.00 

82.00 

82.00 

[Eg]out, kmol/m 1 

0.1 1 10 

0.1102 

0.1091 

0.1101 

T0' r (E a ]o« t , 

kmol/m 1 

9.58 

9.78 

10.38 

10.00 

10 4 [EvW, 

kmol/m' 1 

7.39 

7.50 

7.8 

7.5 

[Edik; ]out. 

kmol/m' 1 

0.1675 

0.1681 

0.1692 

0.1682 


* All the values given above are for DPd - 82.0 and DPf - 40.0 


«tfts *» a-43- 3.91a 




1 Q 4 [E v ] ou t f kmol/frr 1 0 4 [E J ou i kmol/rrr 1 0 4 [E V ] ou t, kmol/m 



Fig. 4.1 [E a | out vs. |EvI „ ut for the feasible chromosomes for different values of 
generation number. N. . The first three generations are shown. NSGA - II used 


• • 








9.4 9.6 9.8 10.0 10.2 10.4 


10 4 [E a ] out ,kmol/m 3 

Fig. 4.3 [E a ] out vs. [Ey] out for the feasible chromosomes for N g = 100. Filled circle - 
Optimal point obtained using NSGA - II. Filled Hexagon- Optimal point obtained using 
NSGA - I. Open circle - Optimal point obtained while minimizing Ii* and I 2 * individually. 
Filled diamond represents optimal point as obtained by Bhaskar et al. (2000b). Filled 
triangle represents the present operating condition of the industrial wiped film reactor. 
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The multiobjective function optimization problem with end-point constraints has 
been solved using NSGA - I also. Figs. 4.4 and 4.5 show all the feasible solutions (i.e., 
those satisfying the end point constraints in Eqs. 4.1b-d) at different values of the 
generation number, N g , when NSGA - I is used. For the same generation number, the 
chromosomes are found less scattered in the [E a ] 0 ut - [E v ] 0 ut space when NSGA - 1 is used. 
A unique optimal solution was obtained for this case, too (for Ng > 50). But NSGA - I 
converges to the unique optimal point in fewer generations than NSGA - II. Also, the 
concentrations of the acid and vinyl end groups in the product at the (unique) optimal point 
were higher than those obtained using NSGA - II. The solution is given by, 



NSGA - 1 

NSGA - II 

[E a ] out , kmol/m 3 

9.78 x 10' 4 

9.58 x 10~ 4 

[E v ]out, kmol/m 3 

7.5077 x \0 A 

7.392 x Iff 4 

DPout 

82.000 

82.000 


The uniqueness of the optimal solution is checked out by performing an optimization using a 
single objective function. The objective functions, T* and h*, were taken, one at a time, and 
optimal solutions were generated under the same conditions as given in Eq. 4.1b-j (using NSGA 
- II). A unique optimal solution was obtained for this case also, but the values of the acid and 
vinyl end group concentrations were not the same as that obtained with the two objective 
function optimization problem. The optimal values of the acid and vinyl end group 
concentrations for the two cases are (using NSGA - II) : 
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Min I,* 

Min h* 

DPoul 

82.00 

82.00 

[EJout, kmol/m' 5 

9.61 x 10 -4 

9.61 x 10- 4 

[E v ]out, kmol/m 3 

7.45 x IQ- 4 

7.45 x IQ- 4 


This is in contrast to the study of Bhaskar et al. (2000b) who obtained a unique optimal solution 
when using a single objective function that is identical to that obtained for a two objective 
function optimization problem. This is quite interesting. It should be emphasized that GA (as 
well as NSGA) is an extremely robust technique [Goldberg (1989)] and should converge to the 
global optimal solution and also give all the optimal solutions where several exist. Clearly, this 
is not happening for either of the two algorithms. 

Fig. 4.3 (N g = 100) and Table 4.2 compare our optimal solutions obtained using NSGA - 
1 and NSGA - II with the values corresponding to the current operation of the industrial wiped 
film reactor and that of Bhaskar et al. (2000b). It is observed that the optimal values of both the 
objective functions obtained with NSGA - II are, indeed, better than the present operating 
conditions of the industrial reactor. Lower values for both the acid and vinyl end group 
concentrations are obtained at a lower temperature. The lower temperature under optimal 
conditions lead to some energy savings. 

The unique optimal solution was found to be affected significantly by the value of the 
seed number, Sr, used, in both NSGA - II and NSGA - 1 code (though unique points were always 
obtained, irrespective of the value of Sr). The results are shown in Fig. 4.6 for the six values of 
Sr given in Table 4.3. It is observed from Fig. 4.6a that both NSGA - II and NSGA - 1 converge 
to diff jrent (thou gh unique) solutions, and both are unable to converge to the best solution in a 
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Table 4.3: Values of S r used for the 5-decision variable problem 


# 

S r 

1 or a* 

0.1298 


0.3219 

2 or b 

(Ref. NSGA - II) 

3 or c 

0.3500 


0.6237 

4 or d 

(Ref. NSGA -I) 

5 or e 

0.6755 

6 or f 

0.9740 


0.8887 [Bhaskar et al 

7 

(2000)] 

8 

Industrial 


* 1 ,2,3,. . . Used for NSGA - II 

a,b,c, ... Used for NSGA - 1 
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single application. The sensitivity of NSGA to one of the computational parameters may be a 
manifestation of premature convergence. This has been alluded to by Goldberg (1989) himself, 
who mentioned that even though GAs are quite efficient in reaching the global-optimum region, 
they are not guaranteed to reach the precise location of the global-optimum point(s), particularly 
for complex problems, and may converge prematurely. 

Another interesting point suggested by Figure 4.6 is that even for near-identical values of 
the two objective functions, as for example, points 1 and 5 (for NSGA - II) corresponding to 
values of S r of 0.1298 and 0.6755, and points a and e (for NSGA - 1) corresponding to values of 
S r of 0.1298 and 0.6755, the decision variables differ significantly. We refer to this as the 
problem of sensitivity of optimal solutions. 

The effect of the weightage factor, W 2 , was studied in detail (using both NSGA - II and 
NSGA - 1). Fig. 4.7 shows the results for six sets of values of Sr (as given in Table 4.3, except 
for Sr = 0.321 9 and Sr = 0.6755 for NSGA - II and Sr = 0.6237 for NSGA - 1), when W 2 is taken 
as 10 4 (instead of 10 8 ). Different (though unique) optimal points are obtained, but the optimal 
value of the acid and vinyl end group concentrations obtained for each of these Sr is lower 
(better) than those obtained with W 2 = 10 8 . Thus, our earlier conjecture that premature 
convergence is associated with the use of inappropriate GA parameters, is confirmed. The 
problem of sensitivity was not observed in this case. The use of a lower value of W 2 , however, 
does not give rise to unique solutions when Sr = 0.3219 and Sr = 0.6755 for NSGA - II and Sr = 
0.6237 for NSGA - I. In these cases, deviations of DP 0U t from DPd (of ± 0.15 for W 2 = 10 4 ) are 
observed. These points are not included in Fig. 4.7. Figs. 4.8, 4.9 and 4.10 show values of [E a ] 0 ut 
vs. [EvU and the corresponding values of the decision variables (P,T,[Sb 203 ], N*, 0*) vs. [E a ] 0 ut 
for these three cases. 
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1U ^aJ out> kmol/m 

Fig. 4.7 Optimum solutions for several values of Sr (Table 4.3) for Eq. 4.1. The value of the 
weightage factor, W 2 , is 10 4 . Point b represents the global optimum solution. Values of the 
decision variables, P and T, for the optimal points are shown in Figs. B, C. (• - NSGA - II, 
o-NSGA-I) 






1 0 4 t E alout. kmol/m 3 

Fig. 4.8 Optimum values of [E a ]out vs. [E v ] 0 ut obtained using NSGA - II for Sr = 0.3219 for 
the five decision variable problem of Eq. 4.1 (Figure a). W 2 = 1 x 10 4 . c-g represent the 
corresponding values of the decision variables (P,T,[Sb 203 ]» N*, 0*). Unique optimal 
solution is not obtained in this case. 



0.044 


e 



9.2 9.4 9.6 9.8 10.0 10.2 10.4 


10 4 [E a ] ou t, kmol/m 3 


Fig. 4.9 Optimum values of [E a ] ou t vs. [E v ] 0 ut obtained using NSGA - II for Sr = 0.6755 for 
the five decision variable problem of Eq. 4.1 (Figure a). W 2 = 1 x 10 4 . c-g represent the 
corresponding values of the decision variables (P,T,[Sb 203 ], N*, 0*). Unique optimal 
solution is not obtained in this case. 
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Fig. 4.10 Optimum values of {E a ]out vs. [E v ] 0 ut obtained using NSGA - 1 for Sr = 0.6237 for 
the five decision variable problem of Eq. 4.1 (Figure a). W 2 = 1 x 10 4 . c-g represent the 
corresponding values of the decision variables (P,T,[Sb 203 ], N*, 0*). Unique optimal 
solution is not obtained in this case. 



We modified the problem in Eq. 4.1 to one involving only four decision variables, and 
kept the temperature at a constant value of 564.0 K. The results obtained using NSGA- 1 and 
NSGA - II (using W 2 = lxlO 4 ) for the six sets of values of Sr (Table 4.3) are shown in Fig. 4.11. 
Different (though unique) local optimal points were obtained for this problem also and these 
points do not form a Pareto set, unlike what was observed by Bhaskar et al. (2000c). Similar to 
what was observed previously, lower values of W 2 , do not give rise to unique solutions when Sr 
= 0.3219 and Sr = 0.3500 for NSGA - II and Sr = 0.3219, 0.3500 and 0.6237 for NSGA - 1. In 
these cases deviations of DP 0U t from DPa (of ±0.15 for W 2 = 10 4 ) are observed. These points are 
not included in Fig. 4.11. Values of [E a ] 0 ut vs. [E v ] 0 ut and the corresponding values of the 
decision variables (P,T,[Sb20 3 ], N*, 0*) vs. [E a ] out for these cases are shown in Figs. 4.12 - 4.16. 

4.4 Summary 

An improved re-tuned model for the finishing reactor in PET manufacture is developed. 
Minimization of the degradation side products, [E a ] ou t and [E v ] 0 ut, was carried out using NSGA — 
I and NSGA - II. P, T, [St^Ch], 0* and N* were used as the decision variables. An end-point 
constraint was imposed on DP out , while the concentration, [E a ] out , was constrained to lie below a 
maximum value. In addition, an allowable range was prescribed for [EdegW A unique solution 
was obtained for this multiobjective optimization problem. No Pareto set of non-inferior 
solutions was obtained. Problems of sensitivity and premature convergence were encountered 
when different values of Sr (randomseed numbers) are used to obtain the optimal solutions. 
When w 2 is reduced from 10 8 to 10 4 , problems of sensitivity and premature convergence were not 
encountered, and optimal values of both the objective functions obtained using NSGA - 1 and 
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10 t E a ] out- kmol/m 

Fig. 4.11 Optimum solutions for several values of Sr (Table 4.3) for the four decision 
variable problem at T = 564.00 K. W 2 = 10 4 . (• - NSGA - II, o - NSGA - 1). 






[Sb 2 0 3 ], wt% 



Fig. 4.11 (Contd...) Values of the decision variables, [Sb 2 0 3 ], 0* and N* for the optimal 
points in Fig. 4.11a for several values of Sr (see Table 4.3) . (• - NSGA - II , o - NSGA - 1). 
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Fig. 4.12 Optimum values of [E a ] 0 ut vs. [E v ] 0 ut obtained using NSGA - II for Sr = 0.3219 
for the four decision variable problem at T = 564.00 K (Figure a), wj = 1 x 10 4 . c-g 
represent the corresponding values of the decision variables (P,T,[Sb20s], N*, 0*). Unique 
optimal solution is not obtained in this case. 





9.2 9.4 9.6 9.8 1 0.0 10.2 10.4 


10 4 [E a ] ou t. kmol/m^ 


Fig. 4.13 Optimum values of [E a ]out vs. [E v ]out obtained using NSGA - II for Sr = 0.3500 
for the four decision variable problem at T = 564.00 K (Figure a). W 2 = 1 x 10 4 . c-g 
represent the corresponding values of the decision variables (P,T,[Sb 203 ], N*, 0*). Unique 
optimal solution is not obtained in this case. 
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Fig. 4.14 Optimum values of [E a l out vs. [E v ] out obtained using NSGA - I for Sr = 0.3219 
for the four decision variable problem at T = 564.00 K (Figure a). W 2 = 1 x 10 4 . c-g 
represent the corresponding values of the decision variables (P,T,[Sb 203 ], N*, 0*). Unique 
optimal solution is not obtained in this case. 
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Fig. 4.15 Optimum values of [E a ] ou t vs. [E v ] out obtained using NSGA - I for Sr = 0.3500 
for the four decision variable problem at T = 564.00 K (Figure a). w 2 = 1 x 10 4 . c-g 
represent the corresponding values of the decision variables (P,T,[Sb 2 03], N*, 0*). Unique 
optimal solution is not obtained in this case. 
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0.044 



10 4 [E a ] ou t, kmol/m 3 



10 4 [E a ] ou t, kmol/m 3 


Fig. 4.16 Optimum values of [E a l out vs. [E v ] out obtained using NSGA - I for Sr = 0.6237 
for the four decision variable problem at T = 564.00 K (Figure a). W 2 = 1 x 10 4 . c-g 
represent the corresponding values of the decision variables (P,T,[Sb 203 ], N*, 0*). Unique 
optimal solution is not obtained in this case. 




NSGA - II are almost identical even with the use of different values of Sr. However, in some 
cases due to deviation of DP 0Ut from DP d , unique optimal points are not obtained. Unique optimal 
points were obtained even when the temperature is kept constant and not used as a decision 
variable (except in some cases, where deviation of DP 0U t from DP d value is observed), in contrast 
to the study of Bhaskar et al. (2000c). 
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Chapter 5 : CONCLUSIONS AND SUGGESTIONS FOR FUTURE 

WORK 

The model of Bhaskar et al. (2000a) has been modified and re-retuned to describe 
industrial data for a wiped film PET reactor. Industrial data involving commercially important 
product properties like DP 0U t and the concentrations of the DEG and acid end groups in the 
product, have been used. The parameters tuned are kLa re f, the Flory-Huggins’ interaction 
parameter, %i, the equilibrium constants, Ki, K 5 and Kg, the frequency factors of three-rate 
constants, k 2 , k 3 and kg, and a, a parameter associated with the effect of the RPM of the agitator 
on ki.a These parameters are tuned using three sets of industrial data. The tuned parameters are 
then used unchanged to predict the concentration of the various side products and the DP ou t for a 
different operating temperature. It is found that the model predicts these values very well. 

Minimization of the degradation side products, [E a ] ou t and [E v ] ou t, was carried out using 
NSGA - I and NSGA - II. P, T, [Sb 2 0 3 ], 0* and N* were used as the decision variables. An end- 
point constraint was imposed on DP 0Ut , while the concentration, [E a ] 0 ut, was constrained to lie 
below a maximum value. In addition, an allowable range was prescribed for [EdegW- Unique 
solutions were obtained for this multiobjective optimization problem using both algorithms. No 
Pareto set of non-inferior solutions was obtained. Problems of sensitivity, premature 
convergence and scatter were encountered when both NSGA - I and NSGA - II were used for 
solving this problem. 

5.1 Future Work 

The following are some of the challenging areas to attempt to understand the PET reactor 
systems: 
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> Kinetic model for PET reactor systems could further be explored to have reliable and 
accurate kinetic data. 

> More work could be done to understand the effect of the Flory-Huggins’ interaction 
parameter (xi) and the mass-transfer parameter (kLa) on the polymer formation of PET. 
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